function [ predprice, vpredprice ] = predictprice_kalman( data )
%PREDICTPRICE_KALMAN Summary of this function goes here
%   Detailed explanation goes here
%%
NN = length(data);

predprice = zeros(1,NN);
predprice(1) = data(1);
predprice(2) = data(2);
predprice(3) = data(3);
vpredprice = predprice;

for j =2:NN
    fprintf('%f\n', j);
    y = data(1:j);
    t = size(y,2);
    Ny = size(y,1);
    Nz = Ny;
    F = eye(Nz);
    H = eye(Ny);
    % observation noise variance is diagonal
    Rt = 1e-5*eye(Ny);
    % transition noise variance
    Vt = .1*eye(Nz);
    a = zeros(Nz,1);
    b = zeros(Ny,1);
    myu = zeros(Nz,1);
    sigma = .1*eye(Nz);
    itermat = [];

    for iter=1:100
       % construct big cov matrix
       var = [Vt zeros(Nz,Ny); zeros(Ny,Nz) Rt];

       % initial values are changed
       z0 = F*myu;
       vz0 = F*sigma*F'+Vt;

       % filtering to get one-step prediction and filtered value
       [logl, pred, vpred, filt, vfilt] = kalcvf(y, 0, a, F, b, H, var, z0, vz0);
       logl = (-2*logl-Ny*log(2*pi))*t;

       % smoothing using one-step prediction values
       [sm, vsm] = kalcvs(y, a, F, b, H, var, pred, vpred);

       % store old parameters and function values
       myu0 = myu;
       F0 = F;
       Vt0 = Vt;
       Rt0 = Rt;
       logl0 = logl;
       itermat = [itermat; iter logl0 F0(:)' myu0'];

       % obtain P*(t) to get P_T_0 and Z_T_0
       % these values are not usually needed
       % See Harvey (1991 p154) or Shumway (1988, p177)
       jt1 = sigma*F'/vpred(:,:,1);
       p_t_0  = sigma+jt1*(vsm(:,:,1)-vpred(:,:,1))*jt1';
       z_t_0  = myu+jt1*(sm(:,1)-pred(:,1));
       p_t1_t = vpred(:,:,t);
       p_t1_t1 = vfilt(:,:,t-1);
       Kt = p_t1_t*H'/(H*p_t1_t*H'+Rt);

       % obtain P_T_TT1. See Shumway (1988, p180)
       p_t_ii1 = (eye(Nz)-Kt*H)*F*p_t1_t1;
       st0 = vsm(:,:,t)+sm(:,t)*sm(:,t)';
       st1 = p_t_ii1+sm(:,t)*sm(:,t-1)';
       st00 = p_t_0+z_t_0*z_t_0';
       cov = (y(:,t)-H*sm(:,t))*(y(:,t)-H*sm(:,t))'+H*vsm(:,:,t)*H';
       for i=t:-1:2
          p_i1_i1 = vfilt(:,:,i-1);
          p_i1_i  = vpred(:,:,i);
          jt1 = p_i1_i1*F'/p_i1_i;
          p_i1_i  = vpred(:,:,i-1);
          if (i>2)
             p_i2_i2 = vfilt(:,:,i-2);
          else
             p_i2_i2 = sigma;
          end
          jt2 = p_i2_i2*F'/p_i1_i;
          p_t_i1i2 = p_i1_i1*jt2'+jt1*(p_t_ii1-F*p_i1_i1)*jt2';
          p_t_ii1 = p_t_i1i2;
          temp = vsm(:,:,i-1);
          sm1 = sm(:,i-1);
          st0 = st0+(temp+sm1*sm1');
          if (i>2)
             st1 = st1+(p_t_ii1+sm1*sm(:,i-2)');
          else
             st1 = st1+(p_t_ii1+sm1*z_t_0');
          end
          st00 = st00+(temp+sm1*sm1');
          cov = cov+H*temp*H'+(y(:,i-1)-H*sm1)*(y(:,i-1)-H*sm1)';
       end

       % M-step: update the parameters
       myu = z_t_0;
       F = st1/st00;
       Vt = (st0-F*st1')/t;
       Rt = cov/t;

       % check convergence
       if norm((myu-myu0)./(myu0+1e-6),inf) < 1e-2 && ...
          norm((F-F0)./(F0+1e-6),inf) < 1e-2 && ...
          norm((Vt-Vt0)./(Vt0+1e-6),inf) < 1e-2 && ...
          norm((Rt-Rt0)./(Rt0+1e-6),inf) < 1e-2 && ...
          abs((logl-logl0)/(logl0+1e-6)) < 1e-3
          break
       end
    end

    
    eval = eig(F0);
    
    eval = [real(eval) imag(eval) abs(eval)];

    var = [Vt zeros(Nz,Ny); zeros(Ny,Nz) Rt];

    % initial values are changed
    z0  = F*myu;
    vz0 = F*sigma*F'+Vt;
    itermat=[];

    % multistep prediction
    [logl, pred, vpred] = kalcvf(y, 1, a, F, b, H, var, z0, vz0);
    predprice(j+1) = pred(end);
    vpredprice(j+1) = vpred(:,:,end);
end

end

